% Algorithm 1

for k=1:100
    for i=1:100
        for j=1:100
           %S1(i,j) = (i*cos(k/100*pi)+j*sin(k/100*pi))/200;
            S1(i,j) = (i+j)/200;
            S2(i,j) = (i*cos(k/100*pi)+j*sin(k/100*pi))/200;
        end
    end
    [Q R] = qr(S1);             % Orthogonal completion of S1
    [U V X C S] = gsvd(S1,S2);  % Generalized SVD
    T = diag(acos(diag(C)));    % Find Theta
    A = U*T*transpose(V);       % Compute velocity matrix
    imagesc(vertcat(horzcat(S1,U,C),horzcat(S2,V,S)));
    colormap(flipud(gray));
    axis equal;
    M1(k) = getframe;
end

movie(M1,1,25);